# Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Mario Orsi
# This file is part of Brahms.
# Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
# Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. 

# This program opens the file "waterPermeationEvents.log" (one of BRAHMS's output files) to compute the water permeability coefficient.
# You might need to issue the command: "chmod a+x computePermeability.py" to make the script executable

# USAGE: cat waterPermeationEvents.log | python computePermeability.py | tee Pw

import sys, string
from math import sqrt
lines = sys.stdin.readlines()

## nWaters = 3655 # number of water molecules in the system
## lateralEdge = 6.228 # box x and y edges [nm]
## normalEdge = 6.461  # box z edge [nm]

nWaters = 4232 # number of water molecules in the DOPC system
#nWaters = 5000 # number of water molecules in the DMPC system
#lateralEdge = 6.2e-3 # [micrometers] DMPC box x and y edges
lateralEdge = 6.81e-3 # [micrometers] DOPC box x and y edges
print "nWaters = %6.0f, lateralEdge = %3.2f nm" % (nWaters, lateralEdge*1e3)
print "***MAKE SURE nWater and lateraledge ARE RIGHT!***"
#normalEdge = 80.   # box z edge [nm]


interfacialArea = lateralEdge**2 
#boxVolume = interfacialArea * normalEdge
#waterConcentration = nWaters / boxVolume NOOOOOOOO!
#waterConcentration = 33.3 # /nm^3
waterConcentration = 33.3e9 # /micrometers^3

crossingWaters = 0
crossingWatersLeftToRight = 0
crossingWatersRightToLeft = 0

brahmsWatPermFile = open( 'waterPermeationEvents.log', 'r' )

for line in lines:
    crossingWaters = crossingWaters + 1
    words = string.split( line )
    #print words[11]
    if float(words[11]) < 0: # if z-coordinate negative
        crossingWatersLeftToRight = crossingWatersLeftToRight + 1
    else:
        crossingWatersRightToLeft = crossingWatersRightToLeft + 1
    if crossingWaters == 1:
        initialTime = string.atof(words[0]);
    finalTime = string.atof(words[0]);
    
brahmsWatPermFile.close

print "total crossingWaters = %d" % crossingWaters
print "crossingWatersLeftToRight = %d" % crossingWatersLeftToRight
print "crossingWatersRightToLeft = %d" % crossingWatersRightToLeft

#timeWindow = finalTime - initialTime
#timeWindow = finalTime
#timeWindow = 1.2e-6 # [s]
timeWindow = 120e-9 # [s]
print "timeWindow = %.2f microseconds" % ( timeWindow *1e6 ) # [microsecond]
print "***MAKE SURE timeWindow IS RIGHT!***"

unidirectionalFlux = crossingWaters / 2. / ( timeWindow )
unidirectionalFluxLeftToRight = crossingWatersLeftToRight  / ( timeWindow )
unidirectionalFluxRightToLeft = crossingWatersRightToLeft  / ( timeWindow )

# compute permeability coefficient [m/s]
P = unidirectionalFlux / ( waterConcentration * interfacialArea )
PLeftToRight = unidirectionalFluxLeftToRight / ( waterConcentration * interfacialArea )
PRightToLeft = unidirectionalFluxRightToLeft / ( waterConcentration * interfacialArea )

#print "The permeability coefficient is %4.2f E-3 cm/s (exp: 7 E-3 cm/s)" % (P * 100000)
print "PLeftToRight = %7.4f micrometers / s " % PLeftToRight
print "PRightToLeft = %7.4f micrometers / s " % PRightToLeft
print "The permeability coefficient is %4.1f micrometers / s (exp: 70 micrometers / s )" % P
Pavg = 0.5 * (PLeftToRight+PRightToLeft)
Psd = sqrt(1./2.*(PLeftToRight**2+PRightToLeft**2)-(Pavg)**2) # std dev
Pse = Psd/sqrt(2) # std error
print "P = %7.4f +- %7.4f micrometers / s " % (Pavg, Pse)
